
# Setup ------------------------------------------------------------------------
  # Options
  options(stringsAsFactors = FALSE)
  # Packages
  library(magrittr)
  library(data.table)
  library(stringr)
  library(lubridate)
  library(lfe)
  library(broom)
  library(stargazer)
  # Set the number of threads for the lfe package
  options(lfe.threads = 4)
  # Directories of joined bill-weather files
  dir_project <- "S:/Edward/Elasticity/"
  #dir_project <- "/Users/edwardarubin/Dropbox/Research/MyProjects/NaturalGas/"
  dir_csv     <- paste0(dir_project, "DataCsv/")
  dir_rds     <- paste0(dir_project, "DataR/")
  dir_figures <- paste0(dir_project, "Figures/")
  dir_output  <- paste0(dir_project, "Output/Tables2WayBillClimate/")
  dir_results <- paste0(dir_rds, "Results/Results20160821/")

# Plain OLS --------------------------------------------------------------------
  # Load the relevant files
  ols_mrg <- readRDS(paste0(dir_results, "ols_mrg.rds"))
  ols_mrg_hdd <- readRDS(paste0(dir_results, "ols_mrg_hdd.rds"))
  ols_base <- readRDS(paste0(dir_results, "ols_base.rds"))
  ols_base_hdd <- readRDS(paste0(dir_results, "ols_base_hdd.rds"))
  # Build results table
  stargazer(ols_mrg, ols_mrg_hdd, ols_base, ols_base_hdd,
    float.env = "table",
    type = "text",
    out = paste0(dir_output, "ols.tex"),
    label = "tab:ols",
    keep = c(1, 3),
    header = F,
    title = "OLS-based elasticity estimates",
    covariate.labels = c("Log(Marginal price)", "Log(Base price)"),
    dep.var.caption = NULL,
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    notes = c(
      "Errors are two-way clustered at the household and year-week.",
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("Zip code's HDD during bill", rep(c("F", "T"), 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
    digits = 4)
  # Clean up
  rm(list = grep("^ols\\_[mb]", ls(), value = T))
  gc()

# First stage: Marginal price; HDD IV ------------------------------------------
  # Load the relevant files
  hdd1_1a <- readRDS(paste0(dir_results, "hdd1_1a.rds")); gc()
  hdd1_1b <- readRDS(paste0(dir_results, "hdd1_1b.rds")); gc()
  hdd1_3a <- readRDS(paste0(dir_results, "hdd1_3a.rds")); gc()
  hdd1_3b <- readRDS(paste0(dir_results, "hdd1_3b.rds")); gc()
  # Calculate first-stage F statistics
  f_stats <- c(
    waldtest(hdd1_1a, ~ east_year_hdd | `east_year_hdd:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(hdd1_1b, ~ east_year_hdd | `thm6_exceed0TRUE` |
      `east_year_hdd:utilitysocal` | `east_year_hdd:thm6_exceed0TRUE` |
      `east_year_hdd:thm6_exceed0TRUE:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(hdd1_3a, ~ east_year_hdd | `east_year_hdd:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(hdd1_3b, ~ east_year_hdd | `thm6_exceed0TRUE` |
      `east_year_hdd:utilitysocal` | `east_year_hdd:thm6_exceed0TRUE` |
      `east_year_hdd:thm6_exceed0TRUE:utilitysocal`,
      type = "cluster")[["F"]]) %>% round(2)
  # Build results table
  stargazer(hdd1_1a, hdd1_1b, hdd1_3a, hdd1_3b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage1_mrg_hdd.tex"),
    label = "tab:stage1 mrg hdd",
    header = F,
    omit = c("ca_year_hdd", "bill_hdd"),
    title = paste("First-stage results:",
      "Instrumenting marginal price with Eastern US heating degree days"),
    covariate.labels = c(
      "Eastern HDD",
      "Simulated inst.",
      "Eastern HDD $\\times$ SoCalGas",
      "Eastern HDD $\\times$ Sim. inst.",
      "Eastern HDD $\\times$ SoCalGas $\\times$ Sim. inst."
      ),
    dep.var.caption = "First-stage results",
    dep.var.labels = "Dependent variable: Log(Marginal price)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("F stat. for instruments", f_stats),
      c("California HDD", rep("T", 4)),
      c("Bill's HDD", rep(c("F", "T"), each = 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
    digits = 4)
  # Clean up
  rm(list = grep("^hdd1_[0-9][ab]", ls(), value = T)); rm(f_stats); gc()

# First stage: Baseline price; HDD IV ------------------------------------------
  hdd1_2a <- readRDS(paste0(dir_results, "hdd1_2a.rds"))
  hdd1_2b <- readRDS(paste0(dir_results, "hdd1_2b.rds"))
  hdd1_4a <- readRDS(paste0(dir_results, "hdd1_4a.rds"))
  hdd1_4b <- readRDS(paste0(dir_results, "hdd1_4b.rds"))
  # Calculate first-stage F statistics
  f_stats <- c(
    waldtest(hdd1_2a, ~ east_year_hdd | `east_year_hdd:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(hdd1_2b, ~ east_year_hdd | `thm6_exceed0TRUE` |
      `east_year_hdd:utilitysocal` | `east_year_hdd:thm6_exceed0TRUE` |
      `east_year_hdd:thm6_exceed0TRUE:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(hdd1_4a, ~ east_year_hdd | `east_year_hdd:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(hdd1_4b, ~ east_year_hdd | `thm6_exceed0TRUE` |
      `east_year_hdd:utilitysocal` | `east_year_hdd:thm6_exceed0TRUE` |
      `east_year_hdd:thm6_exceed0TRUE:utilitysocal`,
      type = "cluster")[["F"]]) %>% round(2)
  # Build results table
  stargazer(hdd1_2a, hdd1_2b, hdd1_4a, hdd1_4b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage1_base_hdd.tex"),
    label = "tab:stage1 base hdd",
    header = F,
    omit = c("ca_year_hdd", "bill_hdd"),
    title = paste("First-stage results:",
      "Instrumenting baseline price with Eastern US heating degree days"),
    covariate.labels = c(
      "Eastern HDD",
      "Simulated inst.",
      "Eastern HDD $\\times$ SoCalGas",
      "Eastern HDD $\\times$ Sim. inst.",
      "Eastern HDD $\\times$ SoCalGas $\\times$ Sim. inst."
      ),
    dep.var.caption = "First-stage results",
    dep.var.labels = "Dependent variable: Log(Baseline price)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("F stat. for instruments", f_stats),
      c("California HDD", rep("T", 4)),
      c("Bill's HDD", rep(c("F", "T"), each = 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
    digits = 4)
  # Clean up
  rm(list = grep("^hdd1_[0-9][ab]", ls(), value = T)); rm(f_stats); gc()

# First stage: Spot price IV (marginal price) ----------------------------------
  # Load relevant files
  spot1_1a <- readRDS(paste0(dir_results, "spot1_1a.rds")); gc()
  spot1_1b <- readRDS(paste0(dir_results, "spot1_1b.rds")); gc()
  spot1_2a <- readRDS(paste0(dir_results, "spot1_2a.rds")); gc()
  spot1_2b <- readRDS(paste0(dir_results, "spot1_2b.rds")); gc()
  spot1_12a <- readRDS(paste0(dir_results, "spot1_12a.rds")); gc()
  spot1_12b <- readRDS(paste0(dir_results, "spot1_12b.rds")); gc()
  # Calculate the first-stage F statistics
  f_stats <- c(
    waldtest(spot1_1a, ~
      spot_price_1week | `spot_price_1week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_1b, ~
      spot_price_1week | `spot_price_1week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_2a, ~
      spot_price_2week | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_2b, ~
      spot_price_2week | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_12a, ~
      spot_price_1week | spot_price_2week |
      `spot_price_1week:utilitysocal` | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_12b, ~
      spot_price_1week | spot_price_2week |
      `spot_price_1week:utilitysocal` | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]]
      ) %>% round(2)
  # Build results table
  stargazer(spot1_1a, spot1_1b, spot1_2a, spot1_2b, spot1_12a, spot1_12b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage1_spot.tex"),
    label = "tab:stage1 spot",
    keep = "spot",
    header = F,
    title = "First-stage results for spot-price based instruments",
    covariate.labels = c(
      "Spot price (1-week lag)",
      "Spot price (1-week lag) $\\times$ SoCalGas",
      "Spot price (2-week lag)",
      "Spot price (2-week lag) $\\times$ SoCalGas"
      ),
    dep.var.caption = "First-stage results",
    dep.var.labels = "Dependent variable: Log(Marginal price)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("F stat. for instruments", f_stats),
      c("Bill's HDD", rep(c("F", "T"), times = 3)),
      c("HH-month FE", rep("T", 6)),
      c("City-year-month FE", rep("T", 6))),
    digits = 4)
  # Clean up
  rm(list = grep("^spot1_[0-9]", ls(), value = T)); rm(f_stats); gc()

# First stage: Spot price IV (baseline price) ----------------------------------
  spot1_1c <- readRDS(paste0(dir_results, "spot1_1c.rds")); gc()
  spot1_1d <- readRDS(paste0(dir_results, "spot1_1d.rds")); gc()
  spot1_2c <- readRDS(paste0(dir_results, "spot1_2c.rds")); gc()
  spot1_2d <- readRDS(paste0(dir_results, "spot1_2d.rds")); gc()
  spot1_12c <- readRDS(paste0(dir_results, "spot1_12c.rds")); gc()
  spot1_12d <- readRDS(paste0(dir_results, "spot1_12d.rds")); gc()
  # Calculate the first-stage F statistics
  f_stats <- c(
    waldtest(spot1_1c, ~
      spot_price_1week | `spot_price_1week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_1d, ~
      spot_price_1week | `spot_price_1week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_2c, ~
      spot_price_2week | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_2d, ~
      spot_price_2week | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_12c, ~
      spot_price_1week | spot_price_2week |
      `spot_price_1week:utilitysocal` | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(spot1_12d, ~
      spot_price_1week | spot_price_2week |
      `spot_price_1week:utilitysocal` | `spot_price_2week:utilitysocal`,
      type = "cluster")[["F"]]
      ) %>% round(2)
  # Build results table
  stargazer(spot1_1c, spot1_1d, spot1_2c, spot1_2d, spot1_12c, spot1_12d,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage1_base_spot.tex"),
    label = "tab:stage1 base spot",
    keep = "spot",
    header = F,
    title = "First-stage results for spot-price based instruments",
    covariate.labels = c(
      "Spot price (1-week lag)",
      "Spot price (1-week lag) $\\times$ SoCalGas",
      "Spot price (2-week lag)",
      "Spot price (2-week lag) $\\times$ SoCalGas"
      ),
    dep.var.caption = "First-stage results",
    dep.var.labels = "Dependent variable: Log(Baseline price)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("F stat. for instruments", f_stats),
      c("Bill's HDD", rep(c("F", "T"), times = 3)),
      c("HH-month FE", rep("T", 6)),
      c("City-year-month FE", rep("T", 6))),
    digits = 4)
  # Clean up
  rm(list = grep("^spot1_[0-9]", ls(), value = T)); rm(f_stats); gc()

# First stage: Simualted spot price IV -----------------------------------------
  # Load relevant files
  sim1_1a <- readRDS(paste0(dir_results, "sim1_1a.rds")); gc()
  sim1_1b <- readRDS(paste0(dir_results, "sim1_1b.rds")); gc()
  # Calculate the first-stage F statistics
  f_stats <- c(
    waldtest(sim1_1a, ~
      spot_price_mrg_sim_iv6 | `spot_price_mrg_sim_iv6:utilitysocal`,
      type = "cluster")[["F"]],
    waldtest(sim1_1b, ~
      spot_price_mrg_sim_iv6 | `spot_price_mrg_sim_iv6:utilitysocal`,
      type = "cluster")[["F"]]
      ) %>% round(2)
  # Build results table
  stargazer(sim1_1a, sim1_1b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage1_spot_sim.tex"),
    label = "tab:stage1 spot sim",
    keep = "spot",
    header = F,
    title = "First-stage results for spot-price based simulated instruments",
    covariate.labels = c(
      "Spot price sim. inst. (1-week lag)",
      "Spot price sim. inst. (1-week lag) $\\times$ SoCalGas"
      ),
    dep.var.caption = "First-stage results",
    dep.var.labels = "Dependent variable: Log(Marginal price)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("F stat. for instruments", f_stats),
      c("Bill's HDD", rep(c("F", "T"), times = 1)),
      c("HH-month FE", rep("T", 2)),
      c("City-year-month FE", rep("T", 2))),
    digits = 4)
  # Clean up
  rm(list = grep("^spot1_[0-9]", ls(), value = T)); rm(f_stats); gc()

# First stage: Lagged days IV --------------------------------------------------
  # Load relevant files
  day1_1a <- readRDS(paste0(dir_results, "day1_1a.rds"))
  day1_1b <- readRDS(paste0(dir_results, "day1_1b.rds"))
  # Get F stat for testing the instrument in each regression (squared t)
  f_stats <- c(
    waldtest(day1_1a, ~
      `log(days_lag1)`,
      type = "cluster")[["F"]],
    waldtest(day1_1b, ~
      `log(days_lag1)`,
      type = "cluster")[["F"]]
      ) %>% round(2)
  # Build results table
  stargazer(day1_1a, day1_1b,
    float.env = "table",
    type = "text",
    out = paste0(dir_output, "stage1_days.tex"),
    label = "tab:stage1 days",
    keep = 1,
    header = F,
    title = "First-stage results for billing days instrument",
    covariate.labels = "Log(No. days in bill)",
    dep.var.caption = "First-stage results",
    dep.var.labels = "Dependent variable: Log(Total bill)",
    style = "qje",
    align = F,
    notes = c("Errors are clustered at the household level.",
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("F stat. for IV", f_stats),
      c("Within-bill HDD", rep(c("F", "T"), each = 1)),
      c("HH-month FE", rep("T", 2)),
      c("City-year-month FE", rep("T", 2))),
    digits = 4)
  # Clean up
  rm(list = grep("^day1", ls(), value = T)); rm(f_stats); gc()

# Second stage: Marginal price; HDD IV -----------------------------------------
  # Load relevant files
  hdd2_1a <- readRDS(paste0(dir_results, "hdd2_1a.rds"))
  hdd2_1b <- readRDS(paste0(dir_results, "hdd2_1b.rds"))
  hdd2_3a <- readRDS(paste0(dir_results, "hdd2_3a.rds"))
  hdd2_3b <- readRDS(paste0(dir_results, "hdd2_3b.rds"))
  # Build results table
  stargazer(hdd2_1a, hdd2_1b, hdd2_3a, hdd2_3b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage2_mrg_hdd.tex"),
    label = "tab:stage2 mrg hdd",
    header = F,
    omit = c("ca_year_hdd"),
    order = c("log", "bill"),
    title = paste("Second-stage results:",
      "Instrumenting marginal price with Eastern US heating degree days"),
    covariate.labels = c(
      "Log(Marginal price) (instrumented)",
      "Within-bill HDD"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("California HDD", rep("T", 4)),
      c("Sim. inst. interactions", rep(c("F", "T"), 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
    digits = 4)
  # Clean up
  # rm(list = grep("^hdd2_[0-9][ab]", ls(), value = T)); rm(f_stats); gc()

# Second stage: Baseline price; HDD IV -----------------------------------------
  hdd2_2a <- readRDS(paste0(dir_results, "hdd2_2a.rds"))
  hdd2_2b <- readRDS(paste0(dir_results, "hdd2_2b.rds"))
  hdd2_4a <- readRDS(paste0(dir_results, "hdd2_4a.rds"))
  hdd2_4b <- readRDS(paste0(dir_results, "hdd2_4b.rds"))
  # Build results table
  stargazer(hdd2_2a, hdd2_2b, hdd2_4a, hdd2_4b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage2_base_hdd.tex"),
    label = "tab:stage2 base hdd",
    header = F,
    omit = c("ca_year_hdd"),
    order = c("log", "bill"),
    title = paste("Second-stage results:",
      "Instrumenting baseline price with Eastern US heating degree days"),
    covariate.labels = c(
      "Log(Baseline price) (instrumented)",
      "Within-bill HDD"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("California HDD", rep("T", 4)),
      c("Sim. inst. interactions", rep(c("F", "T"), 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
    digits = 4)
  # Clean up
  rm(list = grep("^hdd1_[0-9][ab]", ls(), value = T)); rm(f_stats); gc()

# Second stage: Spot price IV (marginal price) ---------------------------------
  # Load relevant files
  spot2_1a <- readRDS(paste0(dir_results, "spot2_1a.rds")); gc()
  spot2_1b <- readRDS(paste0(dir_results, "spot2_1b.rds")); gc()
  spot2_2a <- readRDS(paste0(dir_results, "spot2_2a.rds")); gc()
  spot2_2b <- readRDS(paste0(dir_results, "spot2_2b.rds")); gc()
  spot2_12a <- readRDS(paste0(dir_results, "spot2_12a.rds")); gc()
  spot2_12b <- readRDS(paste0(dir_results, "spot2_12b.rds")); gc()
  # Build results table
  stargazer(spot2_1a, spot2_1b, spot2_2a, spot2_2b, spot2_12a, spot2_12b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage2_spot.tex"),
    label = "tab:stage2 spot",
    order = c(2, 1),
    header = F,
    title = "Second-stage results for spot-price based instruments",
    covariate.labels = c(
      "Log(Marginal price) (instrumented)",
      "Within-bill HDD"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("1-week lagged spot price IV", "T", "T", "F", "F", "T", "T"),
      c("2-week lagged spot price IV", "F", "F", "T", "T", "T", "T"),
      c("HH-month FE", rep("T", 6)),
      c("City-year-month FE", rep("T", 6))),
    digits = 4)
  # Clean up
  rm(list = grep("^spot2_[0-9]", ls(), value = T)); gc()

# TODO: ? Run the 12c and 12d; add to tables
# Second stage: Spot price IV (baseline price) ---------------------------------
  # Load relevant files
  spot2_1c <- readRDS(paste0(dir_results, "spot2_1c.rds")); gc()
  spot2_1d <- readRDS(paste0(dir_results, "spot2_1d.rds")); gc()
  spot2_2c <- readRDS(paste0(dir_results, "spot2_2c.rds")); gc()
  spot2_2d <- readRDS(paste0(dir_results, "spot2_2d.rds")); gc()
  spot2_12c <- readRDS(paste0(dir_results, "spot2_12c.rds")); gc()
  spot2_12d <- readRDS(paste0(dir_results, "spot2_12d.rds")); gc()
  # Build results table
  stargazer(spot2_1c, spot2_1d, spot2_2c, spot2_2d, spot2_12c, spot2_12d,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage2_base_spot.tex"),
    label = "tab:stage2 base spot",
    order = c(2, 1),
    header = F,
    title = "Second-stage results for spot-price based instruments",
    covariate.labels = c(
      "Log(Baseline price) (instrumented)",
      "Within-bill HDD"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("1-week lagged spot price IV", "T", "T", "F", "F", "T", "T"),
      c("2-week lagged spot price IV", "F", "F", "T", "T", "T", "T"),
      c("HH-month FE", rep("T", 6)),
      c("City-year-month FE", rep("T", 6))),
    digits = 4)
  # Clean up
  rm(list = grep("^spot2_[0-9]", ls(), value = T)); gc()

# Second stage: Simualted spot price IV ----------------------------------------
  # Load relevant files
  sim2_1a <- readRDS(paste0(dir_results, "sim2_1a.rds")); gc()
  sim2_1b <- readRDS(paste0(dir_results, "sim2_1b.rds")); gc()
  # Build results table
  stargazer(sim2_1a, sim2_1b,
    float.env = "sidewaystable",
    type = "text",
    out = paste0(dir_output, "stage2_spot_sim.tex"),
    label = "tab:stage2 spot sim",
    order = c("log", "hdd"),
    header = F,
    title = "Second-stage results for spot-price based simulated instruments",
    covariate.labels = c(
      "Log(Marginal price) (instrumented)",
      "Within-bill HDD"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Marginal price)",
    style = "qje",
    align = F,
    notes = c(
      paste("Errors are two-way clustered at the household and",
        "week-city-utility levels."),
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("Bill's HDD", rep(c("F", "T"), times = 1)),
      c("HH-month FE", rep("T", 2)),
      c("City-year-month FE", rep("T", 2))),
    digits = 4)
  # Clean up
  rm(list = grep("^spot1_[0-9]", ls(), value = T)); rm(f_stats); gc()

# Second stage: Lagged days IV -------------------------------------------------
  # Load relevant files
  day2_1a <- readRDS(paste0(dir_results, "day2_1a.rds"))
  day2_1b <- readRDS(paste0(dir_results, "day2_1b.rds"))
  # Build results table
  stargazer(day2_1a, day2_1b,
    float.env = "table",
    type = "text",
    out = paste0(dir_output, "stage2_days.tex"),
    label = "tab:stage2 days",
    order = c(2, 1),
    header = F,
    title = "Second-stage results for billing days instrument",
    covariate.labels = c(
      "Log(Total previous bill) (instrumented)",
      "Within-bill HDD"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    notes = c("Errors are clustered at the household level.",
      "Significance levels: *10\\%, **5\\%, ***1\\%."),
    notes.append = F,
    notes.align = "l",
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("HH-month FE", rep("T", 2)),
      c("City-year-month FE", rep("T", 2))),
    digits = 4)
  # Clean up
  rm(list = grep("^fs[0-9][ab]", ls(), value = T)); rm(f_stats); gc()

# Care second-stage results ----------------------------------------------------
  # Load relevant files
  spot2_1a_care <- readRDS(paste0(dir_results, "spot2_1a_care.rds"))
  spot2_1b_care <- readRDS(paste0(dir_results, "spot2_1b_care.rds"))
  spot2_1c_care <- readRDS(paste0(dir_results, "spot2_1c_care.rds"))
  spot2_1d_care <- readRDS(paste0(dir_results, "spot2_1d_care.rds"))
  # Build results table
  stargazer(spot2_1a_care, spot2_1b_care, spot2_1c_care, spot2_1d_care,
    float.env = "table",
    type = "text",
    out = paste0(dir_output, "stage2_spot_care.tex"),
    label = "tab:stage2 spot care",
    omit = "bill_hdd",
    header = F,
    title = "Elasticity estimates with respect to CARE: Spot price instrument",
    covariate.labels = c(
      "Log(Marginal price)",
      "Log(Marginal price)$\\times$CARE",
      "Log(Baseline price)",
      "Log(Baseline price)$\\times$CARE"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("Within-bill HDDs", rep(c("F", "T"), 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
      omit.table.layout = "n",
    digits = 4)

# Seasons second-stage results -------------------------------------------------
  # Load relevant files
  spot2_1a_summer <- readRDS(paste0(dir_results, "spot2_1a_summer.rds"))
  spot2_1a_winter <- readRDS(paste0(dir_results, "spot2_1a_winter.rds"))
  spot2_1b_summer <- readRDS(paste0(dir_results, "spot2_1b_summer.rds"))
  spot2_1b_winter <- readRDS(paste0(dir_results, "spot2_1b_winter.rds"))
  spot2_1c_summer <- readRDS(paste0(dir_results, "spot2_1c_summer.rds"))
  spot2_1c_winter <- readRDS(paste0(dir_results, "spot2_1c_winter.rds"))
  spot2_1d_summer <- readRDS(paste0(dir_results, "spot2_1d_summer.rds"))
  spot2_1d_winter <- readRDS(paste0(dir_results, "spot2_1d_winter.rds"))
  # Build results table for summer
  stargazer(spot2_1a_summer, spot2_1b_summer, spot2_1c_summer, spot2_1d_summer,
    float.env = "table",
    type = "text",
    out = paste0(dir_output, "stage2_spot_summer.tex"),
    omit = "bill_hdd",
    header = F,
    title = "Summer elasticity estimates: spot price instrument",
    covariate.labels = c(
      "Log(Marginal price)",
      "Log(Baseline price)"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("Within-bill HDDs", rep(c("F", "T"), 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
      omit.table.layout = "n",
    digits = 4)
  # Build results table for summer
  stargazer(spot2_1a_winter, spot2_1b_winter, spot2_1c_winter, spot2_1d_winter,
    float.env = "table",
    type = "text",
    out = paste0(dir_output, "stage2_spot_winter.tex"),
    omit = "bill_hdd",
    header = F,
    title = "Winter elasticity estimates: spot price instrument",
    covariate.labels = c(
      "Log(Marginal price)",
      "Log(Baseline price)"
      ),
    dep.var.caption = "Second-stage results",
    dep.var.labels = "Dependent variable: Log(Consumption, daily avg.)",
    style = "qje",
    align = F,
    omit.stat = c("adj.rsq", "rsq", "f", "ser"),
    add.lines = list(
      c("Within-bill HDDs", rep(c("F", "T"), 2)),
      c("HH-month FE", rep("T", 4)),
      c("City-year-month FE", rep("T", 4))),
      omit.table.layout = "n",
    digits = 4)
